03_valicacion

1. Efecto del tamaño del volumen escaneado

Vamos a estudiar el efecto del volumen escaneado pero ahora en vez de usar el rango como variable que representa ese efecto vamos a evaluar el alto del volumen (perpendicular al suelo).

radial_wind <- ReadNetCDF("../VAD/RMA4/cfrad.20181121_105940.0000_to_20181121_110221.0000_RMA4_0200_02.nc", vars = c("Vda", "azimuth", "elevation")) %>% 
  setDT()

VAD <- with(radial_wind, vad_fit(Vda, azimuth, range, elevation, r2 = 0.8, outlier_threshold = 2.5)) %>% 
  setDT()

VAD %>% 
  .[, spd := sqrt(u^2 + v^2)] %>% 
  .[, `:=`(elev_inf = elevation - 0.5,
         elev_sup = elevation + 0.5,
         range_ini = range - 75,
         range_fin = range + 75)] %>% 
  .[, `:=`(hg_inf = beam_propagation(range_ini, elev_inf)[["ht"]],
           hg_sup = beam_propagation(range_fin, elev_sup)[["ht"]])] %>% 
  .[, dz := hg_sup - hg_inf]

VAD[, height_cut := cut_width(height, 200)] %>% 
  na.omit() %>% 
  ggplot(aes(spd, dz)) +
  geom_point(aes(color = factor(elevation))) +
  geom_smooth(method = lm) +
  scale_color_discrete(name = "Elevación") +
  facet_wrap(~height_cut, scales = "free")

na.omit(VAD) %>% 
  .[, FitLm(spd, dz, se = TRUE), by = .(elevation, cut_width(height, 200))] %>% 
  .[term == "dz"] %>% 
  ggplot(aes(estimate, cut_width)) +
    geom_vline(xintercept = 0, color = "darkgray") +
  geom_point(aes(color = factor(elevation), size = df)) +
  geom_path(aes(group = factor(elevation), color = factor(elevation))) +
  scale_color_discrete(name = "Elevación") +
  scale_size(name = "Grados de \nlibertad") +
  labs(title = "Relación entre la velocidad y la distancia", x = "Pendiente (m/s*km)", y = "Capas (m)") +
  # geom_errorbarh(aes(xmin = estimate - std.error*2, xmax = estimate + std.error*2)) +
  theme(legend.position = "bottom")

## Warning: Removed 3 rows containing missing values (geom_point).

## Warning: Removed 3 rows containing missing values (geom_path).

na.omit(VAD) %>% 
  .[, FitLm(spd, dz, se = TRUE), by = .(cut_width(height, 200))] %>% 
  .[term == "dz"] %>% 
  ggplot(aes(estimate, cut_width)) +
  geom_vline(xintercept = 0, color = "darkgray") +
  geom_point(aes(size = df)) +
  geom_path() +
  scale_size(name = "Grados de \nlibertad") +
  labs(title = "Relación entre la velocidad y la distancia", x = "Pendiente (m/s*km)", y = "Capas (m)") +
  # geom_errorbarh(aes(xmin = estimate - std.error*2, xmax = estimate + std.error*2)) +
  theme(legend.position = "bottom")

## geom_path: Each group consists of only one observation. Do you need to
## adjust the group aesthetic?

## Warning in L$marker$color[idx] <- aes2plotly(data, params, "fill")[idx]:
## number of items to replace is not a multiple of replacement length

ring_fit2 <- function (ring, azimuth, elev, outlier_threshold = Inf) 
{
  nas <- is.na(ring)
  if (sum(nas) == length(ring)) {
    return(list(u = NA_real_, v = NA_real_, r2 = NA_real_, 
                rmse = NA_real_))
  }
  n_outliers <- 1
  while (n_outliers > 0) {
    fit <- stats::lm.fit(cbind(1, cos(azimuth * pi/180), 
                               sin(azimuth * pi/180))[!nas, , drop = FALSE], ring[!nas])
    rmse <- stats::sd(fit$residuals)
    outliers <- abs(fit$residuals) >= outlier_threshold * 
      rmse
    n_outliers <- sum(outliers)
    # browser()   
    ring[!nas][outliers] <- NA
    # ring[!nas][outliers] <-  fit$fitted.values[outliers]
    
    # ring <- rvad:::ring_qc(ring, azimuth, max_na = 0.4)
    nas <- is.na(ring)
    if (sum(nas) == (length(ring) - 3) ) {
      
      return(list(u = NA_real_, v = NA_real_, r2 = NA_real_, 
                  rmse = NA_real_))
    }
    
  }
  r2 <- 1 - stats::var(fit$residuals)/stats::var(ring[!nas])
  coef_cos <- fit$coefficients[2]
  coef_sin <- fit$coefficients[3]
  v <- coef_cos/cos(elev * pi/180)
  u <- coef_sin/cos(elev * pi/180)
  
  # list(azimuth, fit$fitted.values, fit$residuals)
  
  return(list(u = u, v = v, r2 = r2, rmse = rmse))
}

radial_wind[range %in% c(6240, 6390) & elevation == unique(elevation)[4]] %>% 
  ggplot(aes(azimuth, Vda)) +
  geom_line(aes(color = factor(range))) +
  # geom_smooth(aes(color = factor(range)), span = 0.5) +
  labs(title = "VAD para dos anillos particularmente molestos",
       subtitle = "Elevación = 4.57º",
       color = "Rango")

radial_wind[range %in% c(6240, 6390) & elevation == unique(elevation)[4]] %>% 
  .[, ring_fit2(Vda, azimuth, elevation[1], outlier_threshold = 2), by = .(range)] %>% 
  .[, V := metR::Mag(u, v)] %>% 
  .[]

##    range         u         v        r2      rmse        V
## 1:  6240 -3.873376 -11.58958 0.9958829 0.5382776 12.21971
## 2:  6390 -4.105933 -11.36880 0.9958341 0.5287766 12.08753

## Warning: Removed 677 rows containing missing values (geom_path).

## Warning: Removed 695 rows containing missing values (geom_point).

parametros <- CJ(max_na = seq(0.1, 0.3, 0.1),
                 max_consecutive_na = 30,
                 r2_min = seq(0.7, 0.9, 0.05),
                 outlier_threshold = c(Inf, 2, 2.5, 3))